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ABSTRACT 

In a variety of potential flow applications, the modal element method has been shown to signifi- 
cantly reduce the numerical grid, employ a more precise grid termination boundary condition, and give 
theoretical insight to the flow physics. The method employs eigenfunctions to replace the numerical grid 
over significant portions of the flow field. Generally, a numerical grid is employed around obstacles with 
complex geometry while eigenfunctions are applied to regions to regions in the flow field where the 
boundary conditions can easily be satisfied. To handle a wider class of CFD problems, the present paper 
extends the modal element method to include comparison function approximations which do not satisfy 
the governing differential equation. To accomplish this task, a double modal series approximation and 
weighted residual constraints are developed to force the comparison functions to satisfy the governing 
differential equation and to interface properly with the finite element solution. As an example, the method 
is applied to the problem of potential flow in a channel with two-dimensional cylindrical like obstacles. 
The calculated flow fields are in excellent agreement with exact analytical solutions. 


INTRODUCTION 

Analytical-numerical approaches that help reduced computer storage requirement offer signifi- 
cant hardware and operational cost savings. The modal element method is one such analytical -numerical 
strategy to grid reduction. The method employs eigenfunctions to replace the numerical grid over signifi- 
cant portions of the flow field. Generally, a finite element grid is employed around obstacles with com- 
plex geometry while eigenfunctions are applied to regions where the boundary conditions can easily be 
satisfied. At the interface between the solution domains, continuity and gradient constraints couple the 
prescribed analytical and finite element solutions. In electromagnetics, acoustics and potential flow 
problems, the modal element method has been found to significantly reduce the numerical grid and 
provide a more precise grid termination boundary condition. More importantly, the method gives some 
theoretical insight to both the physical problem and the numerical analysis. 

So far, the method has been restricted to linear problems where exact solutions (eigenfunctions) 
to the governing differential equation exist. In acoustics, Astley and Eversman (1981) and Baumeister 
and Kreider (1992) applied the modal element method to duct acoustic and scattering problems while 
Baumeister and Baumeister (1994) applied the method to potential flow over a cylinder. However, for 
more general CFD and acoustic problems, eigenfunction solutions do not exist. Therefore, the modal 
element method requires a wider class of approximate analytical functions to handle the more general 
cases. 

The present paper extends the modal element method to include comparison function approxima- 
tions (Meirovitch, 1967, pg. 140) that satisfy natural boundary conditions but do not satisfy the governing 
differential equation. To accomplish this task, a double modal series approximation is developed. Then, 
weighted residual constraints are developed to force this series to both satisfy the governing differential 
equation in the analytical domain and to interface properly with the finite element domain. 



Emphasis in the paper is on the problem formulation although the analysis focuses on a specific 
geometry. As an example, the method is applied to potential flow in an infinitely long channel with two- 
dimensional cylindrical like obstacles. In the present paper, first the method of analysis and domain 
decomposition are discussed. The development of the analytical double series approximation for the 
potential immediately follows. Next, differential and subdomain interface constraints are developed. The 
finite element solution procedure and the flow geometry are then summarized. Finally, three example 
calculations are presented. 


NOMENCLATURE 

A total dimensionless area of finite element domain 
A< e > area of element e 

A m eigenfunction modal amplitude, equation (32) 

A" m n comparison function modal amplitude, equation (4) 

Um,n comparison function modal amplitude, equation (12) 

b dipole strength, equation (30) 

D solution domain, equation (6) 

e(x,y) error, equation (5) 

G global matrix, equation (26) 

h dimensionless channel height, equation (2) 

M number of elements in finite element domain 
M c number of transverse modes, equation (4) 

M C oef tota * number of modal coefficients used in comparison function, equation (23) 

M pts number of grid points on interface S a used in integration, equation (22) 

m mode number, equation (4) 

m* weight mode number, equation (7) 

N number of nodes in finite element domain 

N c number of exponential terms, equation (4) 

N* e > local linear triangular interpolation function, N fe) (x,y), [Nj^(Xj,y.) = g. . (i = 1,2,3; j = 1,2,3)] 
n exponential mode number, equation (4) 


2 


n 


exponential weight number, equation (7) 

R column vector, right hand side of global matrix equation, equation (26) 

Rj global residual error at node I, equation (24) 

S line interface about finite element domain, see figure 3 

S a entrance plane, figure 3 and equation (24) 

S b exit plane, figure 3 and equation (24) 

S + region exterior to S, equation (13) 

S' region interior to S, equation (13) 

s arc length parameter on S, equation (14) 

U+ dimensionless free stream velocity, equation (2) 

U^f normalizing velocity, equation (2) 

Wj (e) local weight function associated with node I, equation (24) 

W m * n * gl°bal and interface weight function, equation (7) 
x dimensionless axial distance, equation (2) 

x^ starting position of finite element grid, see figure 2 
x Q axial intercept of obstacle, see figure 1 
x out ending position of finite element grid, see figure 2 
y dimensionless transverse distance, equation (2) 

y o height of obstacle, see figure 1 

a axial decay constant, equation (4) 

V 2 Laplace operator 

8jj Kronecker delta (5 jj = 1 for 1 = J; 8 jj = 0 for I ^ J) 

<I> column solution vector, equation (27) 

§ dimensionless potential, equation (2) 
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Subscripts: 

a analytical solution entrance domain, see figure 3 

b analytical solution exit domain, see figure 3 

I global node index in finite element domain 

I s global node index for interface S, equation (17) 

i local element index, discussion of equation (25) 

j local element node index, equation (25) 

Superscript: 

- approximate solution 

# dimensional quantity 

( ) average value 

(e) element value 


METHOD OF ANALYSIS 

The present study adapts modal element comparison functions to compute the potential flow field 
in a duct with two-dimensional cylindrical like obstacles, as shown in figure 1 . The exact shape of the 
obstacle is defined by an infinite row of doublets transverse to a uniform flow (Kirchhoff, 1985). For 
obstacles less than half the height of the duct, its shape is nearly a circular cylinder. For this obstacle an 
exact solution for the flow field exists for validating the theoretical results. The detailed shape of the 
obstacle and the exact solution for the potential will be described later in this report. 

For this configuration, a uniform velocity profile U+ is assumed to exist far upstream. For 
in viscid and irrotational two-dimensional flow, the dimensionless potential equation governs; 

V 2 (|) = 0 (1) 

where the following dimensionless variables are employed: 
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<t> # 

*~U +# h # 

oo 


The superscript # indicates a dimensional variable. The channel height h # is used as the characteristic 
dimension. The superscript + indicates the direction of the velocity in the + x direction. All symbols are 
defined in the nomenclature. 

The common numerical approach to this problem (Chung, 1978 or Hilton and Owen, 1979) is to 
extend the grid far upstream and downstream of the obstacle, as shown in figure 2. Here, the assumption 
of uniform velocity at the grid termination is valid. Generally, a large entrance and exit grid region must 
be selected. The leading edge of the obstacle is defined by x 0 while the inlet and exit positions of the grid 
are defined by x^ and x out respectively. Symmetry about the cylinder was not used allowing the computer 
program to handle a greater variety of problems. 

In contrast to the conventional approach, figure 3 displays the modal element grid for this prob- 
lem. The infinite spatial problem domain is divided into three subdomains, the bounded finite element 
domain, which is characterized by the cylindrical obstacle and the surrounding unbounded uniform 
channel entrance and exit domains. The velocity potential is represented approximately in the grid based 
domain by a finite element solution and will be represented analytically by a comparison trial function 
approximation in the uniform semi-infinite entrance and exit domains. 

The finite element domain contains a nodal grid system that covers the region of complicated 
geometry. Linear triangular elements are used and the subdomain interface is approximated by piecewise 
linear segments. In the finite element domain, an approximate solution for the velocity potential is calcu- 
lated by the Galerkin method. In the analytical domains, which extends to ± infinity, a double series 
expansion for the velocity potential is developed. The modal element method couples the analytical and 
numerical solutions by imposing continuity on the potential and velocity at the interface between the 
subdomains shown in figure 3. This coupling results in a single matrix equation. The series coefficients 
and the potential at the finite element nodes are calculated simultaneously, yielding a global representa- 
tion of the potential field. 


COMPARISON FUNCTION SOLUTION 

Meirovitch (1967) suggest using comparison trial functions (satisfies geometric and natural 
boundary conditions) for approximate solutions to partial differential equations. In general, the trial 
function will depend on the geometry of the problem under consideration. Consider the entrance portion 
of the duct upstream to the duct obstruction as shown in figure 3 for x < x^. For this flow field, the 
present solution will employ comparison trial functions which satisfy the natural boundary conditions 
along the upper wall and symmetry centerline for full cylinder in duct, namely that the transverse gradi- 
ent of the velocity is zero: 
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3<b 

v = — = 0 at y = 0& y = h = l 

dy 


( 3 ) 


Therefore, in the entrance region a of figure 3, the comparison trial function will be assumed in 
the form 


M, 


f N, 


6 = U + x + A~ + 

T a <* o 


ii 


e 

m,n 


n<x(x-x.) 


cos(m k y) 


=\yn=\ 


(4) 


The subscript a indicates that this function applies to spatial domain a of figure 3. The tilde above 
the function indicates that it is an approximate estimate of the solution. The first term represents the 
boundary condition on velocity at - infinity. The second constant term is included since it is a solution 
of the governing equation (1). The terms represented by the double summation will blend the distorted 
potential around the obstacle into the uniform potential upstream. In the summation, the trigonometric 
terms cos (m n y) will be called a mode, similar to the terminology used in acoustics. The modal shapes 
represented by cos (m it y), as seen in figure 4, satisfies the natural boundary condition as given by 
equation (3). 

For generality and ease in later integrations, the x dependance of each transverse mode 
cos(m n y) is assumed exponential in form in equation (4) with the decay factor na. Decay will occur 
since the value of (x - x^) will always be a negative quantity in the exponential term of equation (4), 
since (x - x^) has a negative value in the analytical entrance region shown in figures 1 and 3. For the 
problem under consideration, the transverse modes cos(m n y) will decay in a specific fashion which 
depends on the size of the cylindrical like obstacle shown in figure 1 . The accuracy of equation (4) to 
predict this decay will depend on the total number of chosen exponential terms N c and the chosen 
(known) value of a. 

A total of N c decay exponentials will be available to approximate the decay for each assumed 
transverse mode cos(m n y). In effect, a spectrum of decay factors, as shown in figure 5, is applied to each 
m th mode in equation (4). This spectrum is represented by the summation on n which is associated with 
the exponential terms in equation (4). For an acceptable numerical solution, the analysis must determine 
the proper value of the unknown modal coefficient A n ^ n for each assumed exponential decay factor n a 
shown in figure 5. The numerical constraints to be developed later will force the correct axial decay of 
each mode represented by cos(m it y). 

The double series expansion given by equation (4) does not satisfy the governing differential 
equation (1). Substituting equation (4) into the governing equation (1) yields 


M c N c 

V 2 <i> a = X A m.n[( n2a2-m27t2 ) en a( *~ X “ } cos(m ji y)j = e(x, y) (5) 

m=l n=l 

For the special case with a = it and A^ n = A^ 5 m n where 5 m n is the Kronecker delta symbol, 
equation (4) reduces to the eigenfunction solution and the error e(x,y) is identical to zero. This case will 
be used later to validate the general theory. However, for a general assumed value of a, an error e(x,y) 
exists throughout region a. Consequently, the coefficients in equation (4) must be constrained so that the 
governing equation (1) is satisfied in some approximate manner and the error minimized over the domain. 
The method of weighted residuals is used to constrain the coefficients. The residual error associated with 
equation (4) over the analytical domain a is expressed as 
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R m*,n* = j J W m*,n* V \ d ? <** 

X = -~ y = 0 


m* = 1,2, M c 

n* = 1, 2, N c 


( 6 ) 


where the weight function is defined as 


W * *=e n * a(x * in ^ cos (m*7C y) 

m*,n* 


(7) 


The superscript * designates the particular weight index to differentiate it from the summation indices m 
and n used in equation (4). In contrast to FE weighted residual theory using local weight functions, the 
weight function used here is global in nature acting over the whole x and y domain. 

The residual defined by equation (6) can not be set to zero, because the coefficients in equa- 
tion (4) must also be constrained to match the continuity of potential at the interface S a between the finite 
element and analytical domains. A total residual constaint will be defined later. This total residual will be 
set to zero. In addition, the constant A' 0 need not be constrained since it already satisfies the governing 
differential equation. An interface constraint to be presented later will set the value of A Q . 

Substituting equations (5) and (7) into equation (6) and rearranging yields the following expres- 
sion for the residual 


N c *i„ M c 1 

R m*n* = X Je (n * +n)a(x-Xin) ^A" n (n 2 a 2 -m 2 7i 2 ) jcos(m* ny)cos(m rcy)dy dx 

n=l x=-oo m=l y=0 


m* = 1,2, M c (8) 

n* = 1,2 N c 


Noting that 


i j 

| cos(m * n y)cos(m n y) dy = — (m* = m) 

= 0 (m* m) 


y=0 


(9) 


Substituting equation (9) into equation (8) yields 
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R m*.n* = jZ A m*,n( n2(x2 


n=l 



X = -oo 


e (n*+n)a(x-x i0 ) dx 


m* = 1,2, M c (10) 

n* = 1,2, N c 


The m summation in equation (8) has a value only when m equals the chosen weight value m*. 
Consequently, the orthogonality of the cosine has lead to the elimination of the summation on the m 
modes. The final step is to integrate the simple exponential integral to obtain 


R 



1 ^ (n 2 a 2 - m* 2 7c 2 ) 

~ 2 , m *’ n a(n*+n) 

n=l 


m* = 1,2 M c (11) 

n* = 1, 2 N c 

Shortly, this equation will be combined with an interface equation to effectively constrain equation (4) to 
satisfy both the differential equation and join properly with the finite element solution. 

For the exit duct where x > x out , a similar analytical potential is developed in terms of constant 
B+ as defined by 


M„ N 


= U+x + B+ + 

OO O 




-na(x-x ) 


cos(m n y) 


m=l n = l 


( 12 ) 


which gives a similar residual equation. 


INTERFACE CONDITIONS 

At the interface S between the finite element domain and the analytical domains, the local conti- 
nuity of potential can be expressed as 



(13) 


A residual error is associated with matching the analytical potential in equation (4) to the value of 
the potential in the finite element domain along interface S a in figure 3. This residual can be expressed 
numerically by an integral weighting procedure (Baumeister & Baumeister 1994) as follows: 
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( 14 ) 



y x 

fs a W m*,n*M>-<U ds 

y=0 a 


m* = 1, 2 ,M C 

n* = 1,2, ,N C 

where the weight function is given by equation (7). Again, the superscript * designates the particular 
weight index to differentiate it from the trial function indices’ m and n used in equation (4). 
Substituting equation (4) and (7) into equation (14) and rearranging, yields 


y =1 m c 

R m*,n* = -t u ~ x + A o i j s a cos ( m * n y) d y - X 
y =0 m=1 

y=l 

+ j> s cos(m * n y )4> dy (15) 

y=o * 


m,n 


j — ■ 

j> s cos(m*7t y)cos(m 7t y)dy 


m* = 1,2, M c 

n* = 1,2 ,N C 

The exponential terms have been evaluated at x ^ so that the value of the exponential are identi- 
cal to unity. Thus, the exponentials due not appear in equation (15). Because of the harmonic properties 
of the cosine, the first integral in equation (15) is zero. Again for the second integral, the orthogonal 
properties of the cosine, equation (9), yield nonzero integral values only when m is equal to the weight 
m*. Thus, the summation on m can be eliminated from the above equation to give 


XN„ 


y?i 


a . . = 


m*,n* 


0.5A 


m*,n 


, n=l 


) 


k 

y=0 


cos(m * n y)<|> dy 


(16) 


m* = l,2, ,M C 

n* = l,2,....,N c 

In the integral in equation (16), 4> is a function of y prescribed by the finite element nodal points 
on the interface surface S a It suffices to apply a simple quadrature to obtain acceptable results when 
evaluating equation (16). S a is divided into subintervals centered at points (x Is ,y Is ), which correspond to 
finite element boundary nodes introduced later in the paper. The nodes are evenly spaced on the bound- 
ary with the index I s . Once the number of modes M c has been assumed (as defined in equation (4)) the 
number of nodes on S a is set greater or equal to 6M C to adequately resolve all the harmonic terms in 
equation (16). 
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Applying the midpoint rule to the chosen nodes gives 


+ X cos(m*Jiy,J<|> Is Ay t (17) 

m* = 1,2, M c 

n* = 1, 2,....,N C 

where Ay is the distance between nodes except at the two end points. At the ends, Ay is half the distance 
between nodes. 

Equation (17) is the residual error if the analytical potential and the finite element potential do not 
match on the boundary. Shortly, this interface residual will be combined with the differential equation 
residual to get the final constraint equations for the trial function coefficients. But first, the residual error 
for the lowest plane mode will be determined. 

The method of residuals for the lowest plane mode A^ is written as 

y=J 

R o‘= j S|j ~ = 0 (18) 

y=o * 


R 


m*,n* 


0.5A 


m’.n 


The weight used here is unity. The residual error in this case is set to zero. Substituting in the 
value of ^ from equation (4) into equation (18) and rearranging yields 

y=1 N c M c y=l y=l 

R o a =-[ u ~x in +A;] | s> ds-]T n cos(mny)+ <t> ds = 0 (19) 

y— 0 n=l m=l y =0 y=0 


The first integral is identical to 1 . Because of harmonic properties of the cosine, the second 
integral is identical to zero. The third integral is integrated numerically as before in equation (17). Thus, 
the interface constraint on the lowest order coefficient plane mode A ‘ becomes 

m p* 

a 0 - 5 >> =- u :* in (20) 

i,=i 

A similar set of equations can be written for the outlet b region coefficients. 


TOTAL RESIDUAL CONSTRAINTS 

The weight residual approach will be made to constrain the trial function constants in equation (4) 
so that they satisfy both the differential equation and the interface constraint simultaneously. The domain 
and line residuals are summed and set to zero as follows: 
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( 21 ) 


m* = 1,2, ,M C 

n* = 1,2 N c 

There must be one residual equation for every unknown constant ^ n following the summation 
signs in equation (4). Thus, M c X N c residual equations are required. These M c x N c equations, repre- 
sented by equation (21), simultaneously satisfy both the differential equations and the interface condition 
for the higher order modes on the right hand side of the summation signs in equation (4). Substituting 
equations (11) and (17) into equation (21) and rearranging yields 



(n 2 a 2 -m*Vp 
a(n *+n) 


A m*,n " X COS ( m * 11 y I s ) *I,‘ * y l 5 = ° 


I =1 


m* = 1, 2 ,M C (22) 

n* — 1,2, ,N C 

One additional equation, however, is needed for the unknown lowest order mode coefficient A Q . 
But, the differential equation residual associated with A' is identical to zero. Therefore, the total residual 
is simply equation (20). Thus, the total number of auxiliary residual equations is now defined as M coef 
and is equal to 


M coef =1 + M c XN c (23) 

These equations represented by equations (20) and (22) are later combined with the finite element 
equations to form a matrix system that yields values of all the unknowns <f>, values at the FE nodes as 
well as the comparison function constants A' . An additional set of residual equations similar to equa- 
tion (22) is required to determine the B^ n unknowns as defined in equation (12). The B^ n coefficient 
determine the potential in the exit portion of the duct region shown in figure 3. 


FINITE ELEMENT SOLUTION 

The finite element domain, with total area A, is divided into M discrete triangular elements, A (e) , 
e = 1,2,...,M, defined by N comer nodal points (x,,yj), I = 1,2,...,N. The comer nodes for element area 
A< e ) are denoted (xj^y ,<*>), (x 2 (e) ,y 2 (e) ), and (x 3 < e >,y 3 (e) ). The unknown potential 0, is defined at node I. 
To determine the finite element difference equations for the potential at the nodes, application of the weak 
form of the weighted residuals to equation (1) yields 
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(24) 


0 = R, = 



W r (e) ds . 


(I = 1,2, 3 N) 

Equation (24) is the standard weight residual representation of the nodal difference equations for 
the Laplace equation. A detailed derivation of equation (24) can be found in Baumeister and Baumeister 
(1994) or in standard text books (Chung, 1978 or Hinton and Owen 1979). 

To evaluate the integrals in equation (24), it is necessary to represent <|>(x,y) locally. Let N^ e \ 
j = 1,2,3, be the local linear shape functions for linear triangular elements associated with each comer 
node (defined in Hinton and Owen, 1979, pg. 180), so that 

$ (e) (x,y) = ^Nj e) (x,y)<|>( e) = [n ( c) (x. y )]{<)> (e> }• (25) 

j=i 

Also implement the Galerkin method — when the global index I equals the local node index i 
associated with node (Xj^.y^), let equal the local shape function The global shape function 
W, is assumed to be identically zero for any element where node I does not appear (simple pyramid 
weight approximation); thus, the line integrals in equation (24) vanishes unless node I is on the boundary 
S a orS b . 

The nodal difference equations represented by equation (24) are formed by performing the usual 
element by element summation as presented in introductory Finite Element texts. The natural boundary 
conditions are represented by the gradient of the analytical potential in equation (24). These gradient 
terms can be evaluated directly from equations (4) and (12) in terms of the unknown amplitude coeffi- 
cients A" and B + . This insures continuity of the velocity at the interface between the analytical and finite 
element domains. 


GLOBAL MATRIX EQUATION 


The form of the global matrix equation is obtained by combining the finite element difference 
equations described by equation (24) along with the comparison function “difference” equations (20) 
and (22) and their exit equivalent in region b. The global matrix equation has the usual form. 

[G]{0} = {R} (26) 

where 


{®} - {A 0 , Aj j,Aj 2 , A 1,N C ,A 2,1’'” A 2,N C ’ 

A M c ,1’ A M c ,2’”- A M c ,N c ’ 

< t > l ,< t > 2’ < t*3’ < l > 4’ < t ) 5' , ) > 6’ ^N* 


D + D + Tj+ n+ D+ n + 

a o’ b U )15 U b l,N ’ B 2,l’ ’ B 2,N • 

c 1 ’ c 




( 27 ) 
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The {R} column vector contains the free stream velocity terms present in equation (20) and the 
free stream velocity terms that enter equation (24) through the derivatives of <(> a and 4)^. 

The matrix G has the following general form 


A (11) ! <t> (12) i 0 

"o'T^Tb®' 


equations (20) & (22) 

Finite Element Terms - equation (24) 
exit equivalent of equations (20) & (22) 


( 28 ) 


The top and bottom rows in the matrix represent the contribution from the residuals which force the 
comparison function coefficients to satisfy the governing differential equation and interface matching 
condition, equations (20) and (22) and their exit equivalent. The middle row represents the finite element 
contribution from equation (24). 

The submatrix A (11) is a full M coef x M coef matrix composed of the coefficients of the A Q and 
A mn terms in ec l s - {2{)) and (22 ' ) - The submatrix < t )(12) is a M coef X N matrix composed of the 

coefficients of <t> r in equation (20) and (22). Only the nodes on the boundary line S a have nonzero coeffi- 
cients. The submatrices <|) (32) and B (33) are similar in form resulting from applying continuity at the exit 


interface. 

The submatrix A (21) is an N x M coef matrix composed of the coefficients of A" 0 and A' mn from the 
surface integral in equation (24). For each boundary node, there is a full row of terms in A (21) . For 
interior nodes not on the boundary, the surface integral in equation (24) does not contribute to the differ- 
ence equation. Consequently, there is a full row of zeros in A^ 2 ^ at the interior nodes. A similar inter- 
pretation applies to B< 23) . 

<|>( 22 ) is a sparse, highly banded N x N matrix composed of the coefficients of <t>j resulting from 
the solution of equation (24). The global matrix equation, equation (26), is solved by a standard frontal 
solver program. 

Finally, one additional constraint is required to keep the matrix equation (26) nonsingular. Up to 
this point, the difference equations are not completely independent, as discussed in Baumeister and 
Baumeister (1994). The potential must be given a value (grounded) at some nodal location in the finite 
element grid. For the purposes of this paper. 


<(>1 = 0.0 at x = 0.0 y = 1.0 


(29) 


GEOMETRICAL MODEL AND EXACT SOLUTIONS 

The exact shape of the cylindrical like obstacle is defined by an infinite row of doublets trans- 
verse to a uniform flow (Kirchhoff 1985). Assuming an up stream value of U^ = 1 and a channel of unit 
height, the height y as a function of x of the cylindrical like obstacle in figure 1 is given by 

nb 2 sin(Tcy) ( 3( 

2 cosh(rcx)-cos(7ty) 

Equation (30) is solved numerically by the bisection method for the height y of the obstacle as a function 
x and as a function of the assumed strength of the dipole parameter b. The following two values of b are 
used later in the examples: 
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b x c y 0 

0.5642 -0.5084 0.5000 
1.9020 -1.1552 0.9000 

Recall that the extremes of the obstacle x Q and y 0 are defined in figure 1 . The x Q are negative quantities 
since the origin of the (x,y) coordinate system is directly under the center of the cylindrical like obstacle 
to the right of x Q . For dipole strength b equals 0.5642 the obstacle closely approximates a circular cylin- 
der. In general, for obstacles less than half the height of the duct, the shape is nearly a circular cylinder. 

The advantage of using this obstacle is that an exact solution exists for validating the theoretical 
results as given by (Baumeister and Baumeister, 1994) 

7tb 2 sinh(itx) 

<(> = x + — (3 

2 cosh(nx) - cos(rey) 

Equation (31) will be compared graphically to the numerical examples presented later in the paper. 


RESULTS AND COMPARISONS 

To validate the method, three examples are presented for potential flow over a cylinder like 
obstacle described by equation (30). The first two example considers an obstacle that blocks 50 percent of 
the channel while the third example is concerned with an obstacle blocking 90 percent of the channel. For 
the examples to follow, recall that the channel has a dimensionless height of unity and the mean flow 
velocity is unity. The magnitude of the modal coefficients decrease rapidly for the higher order modes; 
thus, the summation parameters are chosen to be M c = 2 and N c = 4 for a comparison to the exact solu- 
tion. Therefore, as defined in equation (23), M coef = 9. Of course, for problems without an exact solution, 
the number of modes must be increased till the results converge. 

Example 1. — Half Channel Obstruction - a = n 

Consider the potential flow over the cylinder with a b value of 0.5642, x 0 = 0.5082 and 
y D = 0.5000. In this example, the finite element region was placed directly over the obstacle as shown in 
figure 3 so that the start of the analytical domain x^ coincides with the beginning of the obstacle at x Q . As 
shown in figure 3, the finite element domain extends from -0.5082 to 0.5082. Thirteen nodes were 
used on the interface and eleven nodes along the surface of the obstacle for a total of 94 nodes and 142 
elements, as shown in figure 6. 

In figure 6, the velocity potential is plotted as a function of x along the upper wall, y = 1 . 

The dashed line represents the potential on the duct wall without an obstacle. The modal element solu- 
tions (hollow boxes) are compared to the exact solution (solid line) given by equation (31). For 
—0.5082 < x < 0.5082, the values of potential at the finite element nodes are used to generate the solution. 
Also in figure 7, the numerical solution is generated from equations (4) and (12) using the numerically 
determined modal coefficients. Six separate values of the potential were calculated in each analytical 
region as shown in figure 7. Clearly, the modal element method gives good agreement with the exact 
results. Figure 8 shows the contour plot of the potential in both the analytical and finite element regions. 
The dash line in figure 8 shows the streamlines. Again, good agreement exists between the exact and the 
modal element results. 

With a = n in this first example, the comparison function solution, equation (4), has imbedded the 
exact eigenfunction solution which has been shown to be (Baumeister and Baumeister, 1994): 
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From equation (31) the exact values of the modal coefficients in equations (4) and (32) have been 
determined to have the following values: 


A = 


O 



(33) 


A = A 8 with A = -Jib 2 e m,n,n (34) 

mn m mn ra 

For this symmetrical obstacle the B+ modes are identical in magnitude to the A^ modes but of 
opposite sign. Because of the sharp decay predicted by equation (34), only the first two higher order 
modes need be calculated to obtain accurate results for the example problems considered. 

To match the exact solution with a = 7t in this first example, the method of weighted residuals 
should set the coefficients to zero when m * n. This is indeed the case. The numerically calculated modal 
coefficients and the exact coefficients are: 


Numerical Solution, equation (4) Exact Solution, equations (33) & (34) 


Aq = -0.487 
Aj j = -0.196 
A, 2 = +0.603x1 0 -4 
Aj 3 = -0.867X10- 4 
Aj 4 = +0.394x1 O' 4 
A 21 = +0.246x1 0“ 5 
A *2 = -0.402x1 0 _1 
A '2 3 = -0.186x10^ 
A 2 ' 4 = -0.858x1 0" 5 


Aq = -0.50002 
Aj , = -0.20247 
Aj’ 2 = +0.0 
A | ^ = +0.0 
A" m = +0.0 

A 2 j — +0.0 

A 22 = -0.40994x1 0 _1 
A 2 ’ 3 = +0.0 
A 2 ’ 4 = +0.0 


As seen above, the numerically calculated and exact modal coefficients are in good agreement. This 
agreement accounts for the excellent match between exact theory (solid line) and the numerical results 
(square symbols) in the analytical regions shown in figures 7 and 8. The agreement occurs in both the 
analytical entrance and exit regions as well as in the finite element region direct on top of the obstacle. 


Example 2. — Half Channel Obstruction - a = 2.5 

Example 2 has the same geometry and parameters as example 1 except a = 2.5. In contrast to 
example 1, for ct = 2.5 the exact eigenvalue solution is not imbedded in the comparison function solution. 
In this case, only an approximate solution is possible. However, the numerical results are still in very 
good agreement with the exact results. To the eye, plots of the numerical results appear exactly as shown 
in figure 7 for the potential along the upper wall and in figure 8 for the contours of the potential. Thus, 
these numerical results are not presented in figure form. An explanation for this excellent agreement can 
be found by examining the coefficients of the comparison function expansion, equation (4). 

The numerically calculated modal coefficients and the exact coefficient are: 


15 



Numerical Solution, equation (4) 

Aj, = -0.484 
Aj , = -0.9489x10-* 

Aj 2 = -0.2041 
Aj' 3 = +0.1553 
Aj 4 = -0.5326x10"* 

Aj’, = +0.7332x1 0" 3 
A' 22 = -0.1516x10-* 

A 23 = -0-3056x10-* 

Aj 4 = +0.5649x1 0" 2 

The numerical value of lowest order mode Aj, is good agreement with the theoretical value. For 
the first transverse mode (m =1), the product of the modal coefficient and its exponential damping factor 
are plotted in figure 9 for various axial mode numbers n. The abscissa is the distance from the interface 
between the finite element region and the analytical inlet region a. In addition, the solid line without 
symbols is a plot of the value of the axial component of the exact eigenfunction given in equation (32). 

An inspection of figure 9 reveals that the sum of the n axial components approximates the value of the 
axial component of the eigenfunction. This is shown explicitly in figure 10. The solid numerical points 
represent the sum of the curves shown in figure 9. Thus, even with a very coarse series approximation, the 
weighted residual approach adjusts the coefficients of the comparison function to match the exact 
eigenfunction. 

Example 3. — Large Channel Obstruction - a = 2.5 

Consider the potential flow over the cylinder with a b value of 1 .9020, x Q = -1 . 1 552 and 
y G = 0.9000. In this case, 90 percent of the channel has an obstruction as illustrated in figure 1 1 . The grid 
was extended slightly in front of the obstacle to better resolve the steep slope of the obstacle near x c . In 
this case, x^ = -1.3984 and x out = 1 .3984. Twenty one nodes were used on the interface and 132 nodes 
along the upper channel wall for a total of 1272 nodes and 2212 elements, as shown in figure 11. 

In figure 12, the velocity potential is again plotted as a function of x along the upper wall, y = 1. 
The dashed line again represents the potential on the duct wall without an obstacle. The modal element 
solutions (hollow boxes) are compared to the exact solution (solid line) given by equation (31). In fig- 
ure 12, for -1.3984 < x <1.3984, the values of potential at the finite element nodes are used to generate 
the solution. Also in figure 12, the numerical solution is generated from equation (4) using the numeri- 
cally determined modal coefficients. Similar to figure 8, the contour plots of the numerical potentials 
were in very good agreement with the theoretical values. 

The numerically calculated modal coefficients and the exact coefficients are as follows: 

Numerical Solution, equation (4) Exact Solution, equation (33) 

A 0 =-0.5671x10* Aj> = -0.5683x10* 

Aj , =-0.6801x10"* 

Aj’ 2 =-0.1465 
Aj’ 3 =+0.1115 
Aj 4 =-0.3828x10"* 

Aj, =+0.3633x10^ 

A- 22 = -0.6749X10" 3 
A - 2 ’ 3 =-0.1 250xl0~ 2 
A- 24 =+0.2165xl0- 3 


Exact Solution, equations (33) 
Aj, = -0.50002x10 
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CONCLUDING REMARKS 


The modal element method using comparison functions has been presented. The method has been 
applied to potential flow over a two-dimensional cylindrical like obstacle. In the example problems, the 
total flow domain is broken into three subdomains that are patched together. The potential field is repre- 
sented by a finite element solution in the irregular subdomain next to the obstacle and by a comparison 
function expansion in the unbounded entrance and exit ducts. By the method of weighted residuals, the 
comparison function is simultaneously constrained to satisfy both the governing Laplace equation and the 
continuity of potential and velocity across the interface between the subdomains. The combined numeri- 
cal and analytical results show excellent agreement with the corresponding exact solutions. The method 
is robust in that relatively coarse series approximations have developed very accurate results. 
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Figure 2. — Conventional finite element discretization for flow around cylindrical object 
in channel. 
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Figure 3. — Modal element finite element discretization for flow around cylindrical 
object in channel. 
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3 dal element finite element discretization for flow around 
object in channel with y 0 = 0.5. 
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X axial coordinate 


Figure 7. — Effect of a half channel obstruction on the potential 
along the upper wall (a = tt, M c = 2, N c = 4 with 94 nodes 

and 1 42 elements). 
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Figure 8. — Contour plots of the potential in both the finite element and the analytical regions for half the channel 
obstruction (a = tt, M c = 2, N c = 4 with 94 nodes and 142 elements). 
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Figure 9. — Axial components of transverse mode m = 1 for 50 percent channel blockage, alpha - 2.5. 


21 



22 


Figure 1 1 . — Modal element finite element discretization for flow around large cylindrical object in channel, 



x axial coordinate 


Figure 1 2. — Effect of a large channel obstruction on the potential 
along the upper wall (1 272 nodes and 221 2 elements). 
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